source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

# AID --------
## SLE ------------
### Perez -------------
sobj <- read_zstd_rds('~/append-ssd/alaria2/GSE174188.perez.SLE.zst.rds')

sobj <- sobj |>
  filter(cov %in% c('B cells','Plasmablasts'))

sobj <- sobj |>
  mutate(orig.ident = ind_cov_batch_cov, ct_cov = NULL) |>
  quick_process_seurat(skip_norm = T, leiden = F)

sobj <- sobj |>
  FindClusters(resolution = .5)

sobj |>
  DimPlot(cols = DiscretePalette(36), label = T)

sobj |>
  DimPlot(group.by = 'cov')

b_cell_markers <- list(
  Naive_B = c("MS4A1", "CD19", "TCL1A", "FCER2", "CCR7", "IGHD", "IGHM"),
  Transitional_B = c("CD24", "CD38", "CD10", "IGLL1", "VPREB1"),
  Memory_B = c("CD27", "TNFRSF13B", "CR2"),
  Plasmablast = c("MKI67"),
  Plasma_Cell = c("SDC1", "XBP1", "PRDM1", "MZB1", "JCHAIN", "IGHG1", "IGHA1"),
  Double_Negative_B = c("FCRL5", "ITGAX", "TBX21", "FCGR3B")
)

sobj |>
  DotPlot(b_cell_markers, cols = 'RdBu', cluster.idents = T) +
  RotatedAxis()

sobj |>
  DotPlot(cc.genes.updated.2019, cols = 'RdBu', cluster.idents = T)

sobj <- sobj |>
  mutate(manual_fine = case_when(
    seurat_clusters == 7 ~ 'DN B cells',
    seurat_clusters == 9 ~ 'Plasma cells',
    seurat_clusters == 11 ~ 'Plasmablasts',
    seurat_clusters %in% c(5,1) ~ 'Memory B cells',
    seurat_clusters == 3 ~ 'Transitional B cells',
    .default = 'Naive B cells'
  ) |> fct_relevel('Transitional B cells', 'Naive B cells',
                   'Memory B cells', 'Plasmablasts',
                   'Plasma cells', 'DN B cells'))

sobj |>
  DimPlot(group.by = 'manual_fine', cols = 'Paired') +
  ggtitle('B cell subsets')

sobj <- sobj |>
  mutate(Status = fct_relevel(Status, 'Healthy', 'Managed', 'Treated'))

sobj |>
  write_zstd_rds('mission/sos1/perez_SLE_B_cell.rds')

#### SLE vs HC -------------
sobj |>
  DotPlot2d('SOS1', SLE_status, manual_fine) +
  labs(x = 'Group', y = 'Cell type', subtitle = 'GSE174188',
       title = 'SOS1 in SLE PBMC B cell subsets')

g1 <- last_plot()

g1 |>
  pluck('data') |>
  mutate(group.y = fct_relevel(group.y, 'Transitional B cells', 'Naive B cells',
                               'Memory B cells', 'Plasmablasts',
                               'Plasma cells', 'DN B cells') |>
           fct_rev()) |>
  BubblePlot(d2 = T, size = c(1,6)) +
  labs(x = 'Group', y = 'Cell type', subtitle = 'GSE174188',
       title = 'SOS1 in SLE PBMC B cell subsets')

sobj |>
  FindMarkersAcrossVar(split.by = 'manual_fine', group.by = 'SLE_status',
                       ident.1 = 'SLE', features = 'SOS1')

sos1_bysample <- sobj |>
  DotPlot2d('SOS1', orig.ident, manual_fine) |>
  pluck('data')

sos1_bysample <- sobj |>
  distinct(orig.ident, Status) |>
  mutate(group.x = as.character(orig.ident), .keep = 'unused') |>
  inner_join(sos1_bysample) |>
  as_tibble()

sos1_bysample |>
  mutate(sle = ifelse(Status == 'Healthy', 'Healthy', 'SLE')) |>
  filter(str_detect(group.y, 'Memo')) |>
  ggplot(aes(sle, avg.exp, fill = sle)) +
  stat_mean(geom = 'col') +
  ggbeeswarm::geom_beeswarm(cex = 3, corral.width = .5, corral = 'wrap',
                            size = 1) +
  facet_wrap(~group.y, scales = 'free_y') +
  stat_compare_means(method = 't.test',
                     comparisons = list(c('Healthy', 'SLE'))) +
  labs(x = 'Group', y = 'Average expression', fill = 'Group',
       title = 'SOS1 expression in SLE',
       subtitle = 'GSE174188') +
  theme_pubr() +
  scale_y_continuous(expand = expansion(mult = c(.05, .1))) +
  scale_fill_hue(direction = -1)

#### Flare vs HC -----------
sobj |>
  DotPlot2d('SOS1', Status, manual_fine) +
  labs(x = 'Group', y = 'Cell type', subtitle = 'GSE174188',
       title = 'SOS1 in SLE PBMC B cell subsets')

g2 <- last_plot()

g2 |>
  pluck('data') |>
  mutate(group.y = fct_relevel(group.y, 'Transitional B cells', 'Naive B cells',
                               'Memory B cells', 'Plasmablasts',
                               'Plasma cells', 'DN B cells') |>
           fct_rev(),
         group.x = fct_relevel(group.x, 'Healthy', 'Managed', 'Treated')) |>
  BubblePlot(d2 = T, size = c(1,6)) +
  labs(x = 'Group', y = 'Cell type', subtitle = 'GSE174188',
       title = 'SOS1 in SLE PBMC B cell subsets')

sos1_bysample |>
  filter(str_detect(group.y, 'Plasma')) |>
  ggplot(aes(Status, avg.exp, fill = Status)) +
  stat_mean(geom = 'col') +
  ggbeeswarm::geom_beeswarm(cex = 3, corral.width = .5, corral = 'wrap',
                            size = 1) +
  facet_wrap(~group.y, scales = 'free_y') +
  stat_compare_means(method = 't.test',
                     comparisons = list(c('Healthy', 'Managed'),
                                        c('Healthy', 'Treated'),
                                        c('Healthy', 'Flare'))) +
  labs(x = 'Group', y = 'Average expression', fill = 'Group',
       title = 'SOS1 expression in SLE',
       subtitle = 'GSE174188') +
  theme_pubr() +
  scale_y_continuous(expand = expansion(mult = c(.05, .15))) +
  scale_fill_hue(direction = -1)

### small one ----------------
sobj <- read_rds('DE_cells/data/sle_sobj.rds')

sobj |>
  DotPlot2d('SOS1', orig.ident, monaco_label) +
  labs(x = 'Sample', y = 'Cell type', subtitle = 'GSE162577 PBMC')

## MS ------------
sobj <- read_rds('DE_cells/data/ms_sobj.rds')

sobj |>
  DotPlot(pbmc_markers, group.by = 'monaco_label', cols = 'RdBu')

sobj <- sobj |>
  mutate(group = ifelse(str_starts(orig.ident, 'MS'), 'MS', 'HC'),
         tissue = ifelse(str_ends(orig.ident, 'CSF'), 'CSF', 'PBMC'))

sobj |>
  filter(tissue == 'PBMC') |>
  DotPlot2d('SOS1', group, monaco_label) +
  labs(x = 'Sample', y = 'Cell type', title = 'SOS1 in MS PBMC',
       subtitle = 'GSE138266')

sobj |>
  filter(tissue == 'CSF') |>
  FindMarkersAcrossVar(split.by = 'monaco_label', group.by = 'group',
                       ident.1 = 'MS', features = 'SOS1')

sobj |>
  filter(tissue == 'CSF') |>
  DotPlot2d('SOS1', group, monaco_label) +
  labs(x = 'Sample', y = 'Cell type', title = 'SOS1 in MS CSF',
       subtitle = 'GSE138266')

ms_sos1 <- sobj |>
  DotPlot2d('SOS1', orig.ident, monaco_label) |>
  pluck('data')

ms_sos1 |>
  mutate(group = ifelse(str_starts(group.x, 'MS'), 'MS', 'HC'),
         tissue = ifelse(str_ends(group.x, 'CSF'), 'CSF', 'PBMC'),
         group.y = ifelse(str_detect(group.y, 'Switched'), 'Memory B cells', group.y)) |>
  filter(str_detect(group.y, 'B')) |>
  ggplot(aes(group, avg.exp, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .2) +
  facet_wrap(vars(tissue, group.y), scales = 'free_y', nrow = 1) +
  theme_pubr() +
  scale_fill_hue(direction = -1) +
  labs(title = 'SOS1 expression in MS B cells', y = 'Average expression',
       subtitle = 'GSE138266')

# pSS ------------
sobj <- read_rds('DE_cells/data/pss_sobj.rds')

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'pSS|HC'))

sobj |>
  DotPlot2d('SOS1', group, monaco_label) +
  labs(x = 'Group', y = 'Cell type', title = 'SOS1 in pSS PBMC',
       subtitle = 'GSE157278')

sobj |>
  FindMarkersAcrossVar(split.by = 'monaco_label', group.by = 'group',
                       ident.1 = 'pSS', features = 'SOS1')

sobj |>
  DotPlot2d('SOS1', monaco_label, orig.ident)

pss_sos1 <- last_plot() |>
  pluck('data')

pss_sos1 |>
  mutate(group = str_extract(group.y, 'pSS|HC')) |>
  filter(str_detect(group.x, 'Plasmab|B')) |>
  ggplot(aes(group, avg.exp, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .2) +
  facet_wrap(~group.x) +
  theme_pubr() +
  scale_fill_hue(direction = -1) +
  labs(title = 'SOS1 expression in pSS B cells', y = 'Average expression',
       subtitle = 'GSE157278')

## RA -----------
sobj <- read_rds('/home/supervisor/mist/NaviLau/RA_zhangxuan/RA_HC_merged.rds')

sobj <- JoinLayers(sobj)

sobj <- sobj |>
  mutate(barcode = str_extract(.cell, '[A-Z]+'),
         unique_bc = str_c(barcode, '_', sample))

colnames(sobj) <- sobj$unique_bc

sobj <- sobj |>
  mutate(orig.ident = sample, barcode = NULL, unique_bc = NULL) |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

Idents(sobj) <- sobj$orig.ident

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj |> VlnPlot('nCount_RNA', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10, nCount_RNA < 40000)

sobj <- sobj |>
  NormalizeData()

sobj |> as_tibble() |>
  write_source_csv('ra.pbmc.metadata')

smeta <- read_csv('mission/SLE_TRPM2_MfMo/results/ra.pbmc.metadata.csv')

sobj <- sobj |>
  left_join(smeta)

sobj |> DotPlot(c(pbmc_markers,'SOS1'), group.by = 'manual_main',
                cols = 'RdBu', cluster.idents = T)

sobj <- sobj |>
  filter(!is.na(manual_main))

sobj <- sobj |>
  mutate(group = genotype, genotype = NULL)

sobj |>
  write_zstd_rds('~/append-ssd/geo_array/RA/HRA000155_RA.zst.rds')

sobj |>
  DotPlot2d('SOS1', group, manual_main) +
  labs(x = 'Group', y = 'Cell type', title = 'SOS1 in RA PBMC',
       subtitle = 'HRA000155')

sobj |>
  FindMarkersAcrossVar(split.by = 'manual_main', group.by = 'group',
                       ident.1 = 'RA', features = 'SOS1')

ra_sos1 <- sobj |>
  DotPlot2d('SOS1', orig.ident, manual_main) |>
  pluck('data')

ra_sos1 |>
  mutate(group = str_extract(group.x, 'HC|RA')) |>
  filter(str_detect(group.y, 'B|Plas')) |>
  ggplot(aes(group, avg.exp, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .2) +
  facet_wrap(~group.y) +
  theme_pubr() +
  scale_fill_hue(direction = -1) +
  labs(title = 'SOS1 expression in RA B cells', y = 'Average expression',
       subtitle = 'HRA000155')
